%% task6_2_3_1.m

clear;
ConstantDefinition();
load("..\data\wave2fft.mat", "overtone");
load("..\data\ConstantDefinition.mat", "SampRate");

dfh_score = [3, 5, 1;
    3, 5, 0.5;
    3, 6, 0.5;
    3, 2, 2;
    3, 1, 1;
    3, 1, 0.5;
    2, 6, 0.5;
    3, 2, 2];

overtone = overtone ./ overtone(1);
global ove;
ove = overtone;

mus = SoundSynthesis6_2_3(dfh_score, 120);
sound(mus, SampRate);
% res = SingleFre6_2_3(440, SampRate, 15, 2);
% plot(res);

%% ConstantDefinition.m

function ConstantDefinition()
% ConstantDefinition 设置全局常量并保存在ConstantDefinition.mat里，用到时load即可

    % F调唱名对应频率
    % 行：大字组，小字组，小字一组，小字二组
    FRE_F = [
             87.31, 98, 110, 116.54, 130.81, 146.83, 164.81;
             174.61, 196, 220, 233.08, 261.63, 293.66, 329.63;
             349.23, 392, 440, 466.16, 523.25, 587.33, 659.25;
             698.46, 783.99, 880, 932.33, 1046.50, 1174.66, 1318.51;
             ];

    % used in function SingleSound
    PHI = 0;                % 初始相位
    TurnPoint1 = 0.25;      % 转折点
    TurnPoint2 = 0.35;
    TurnPoint3 = 0.7;
    Arate = 1.5;            % 最大值占平稳时幅值的比率

    HarmAmp1 = 0.2;         % 第一泛音列幅值占基音的幅值比例
    HarmAmp2 = 0.3;         % 第二泛音列幅值占基音的幅值比例

    % used in function SoundSynthesis
    SampRate = 44100;       % 采样率
    OverlapRate = 0.4;      % 重叠率，与前一个音有所重叠

    save("..\data\ConstantDefinition.mat", "FRE_F", "PHI", "TurnPoint1", "TurnPoint2", "TurnPoint3", "Arate", ...
        "SampRate", "OverlapRate", "HarmAmp1", "HarmAmp2");
end

%% function SoundSynthesis
function res = SoundSynthesis6_2_3(musicScore, bpm)
% SoundSynthesis 将SingleSound得到的许多单音按一定的重叠比率整合为一条音频
%   musicScore: 乐谱参数，默认F调，以下方的fre_f为准。
%               size: [n 3] n行，1,2列分别为小字几组和对应唱名，3行为所占拍数
%               （后期可扩展，实现每个音单独设置其相应参数）（ps：什么UTAU）
%   OverlapRate: 重叠率，后一个音与前一个音的衰减部分比率
%   bpm: 每分钟拍数，衡量音乐速率
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[:,1]

    dur = 60 / bpm;     % 每拍所占时长（单位：秒）
    load("..\data\ConstantDefinition.mat", "FRE_F", "TurnPoint3", "SampRate", "OverlapRate");
    offset = 0;         % 重叠时向前偏移值
    offset_t = 0;       % 重叠时间（单位：秒）
    i = 1;              % 循环变量
    res = 0;
    for index = musicScore'
        sou = SingleSound6_2_3(FRE_F(index(1),index(2)), SampRate, 15, ...
                index(3) * dur + offset_t);
        i = i - offset;
        j = 1;
        while i <= size(res,1) && i > 1
            res(i, 1) = res(i) + sou(j);
            j = j + 1;
            i = i + 1;
        end
        while j <= size(sou,1)
            res(i, 1) = sou(j);
            j = j + 1;
            i = i + 1;
        end
        offset = ceil(OverlapRate * (1 - TurnPoint3) * size(sou, 1));
        offset_t = OverlapRate * (1 - TurnPoint3) * (index(3) * dur + offset_t);
    end
    res(i, 1) = 0;
end

%% SingleSound
function res = SingleSound6_2_3(fre, SampRate, amplitude, durat)
% SingleSound: 单个乐音生成器（叠加泛音列）
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    global ove;
    res = SingleFre6_2_3(fre, SampRate, amplitude, durat);
    for i = 2 : length(ove)
        res = res + ove(i) * SingleFre6_2_3(fre*i, SampRate, amplitude, durat);
    end
end

function res = SingleFre6_2_3(fre, SampRate, amplitude, durat)
% SingleFre: 单频率乐音生成器
% fre: 乐音频率
% SampRate: 采样率
% amplitude: 乐音幅值（持续时幅值）
% durat: 乐音时长（单位：秒）
% res: 返回以正弦波为基础的幅值发生衰减的采样后的离散值
%      大小[时长*采样率 + 1， 1]

    res = zeros(ceil(durat * SampRate), 1);
    load("..\data\ConstantDefinition.mat", "PHI", "TurnPoint1", "TurnPoint2", "TurnPoint3", "Arate");
    ind = 1;            % 循环变量
    for t = 0:1/SampRate:durat
        if t <= TurnPoint1 * durat
            % res(ind,1) = (Arate * amplitude / TurnPoint1 / durat * t) * sin(fre * t + PHI);
            res(ind,1) = (1-10^(-t))/(1-10^(-TurnPoint1 * durat)) * Arate * amplitude * sin(fre * t + PHI);
        elseif t <= TurnPoint2 * durat
            res(ind,1) = (amplitude*(Arate-1)*(exp(-(t-TurnPoint2*durat))-1)/(exp(-(TurnPoint1*durat-TurnPoint2*durat))-1)+amplitude)*...
                sin(fre * t + PHI);
            % res(ind,1) = (-(Arate - 1) * amplitude / (TurnPoint2 - TurnPoint1) / durat * t + ...
            %     (Arate * TurnPoint2 - TurnPoint1) * amplitude / (TurnPoint2 - TurnPoint1)) * ...
            %     sin(fre * t + PHI);
        elseif t <= TurnPoint3 * durat
            res(ind,1) = amplitude * sin(fre * t + PHI);
        else
            res(ind,1) = amplitude*(10.^((1/Arate)-(t-TurnPoint3 * durat)/(durat-TurnPoint3 * durat)*(1/Arate))-1)/(10.^(1/Arate)-1) * sin(fre * t + PHI);
            % res(ind,1) = (amplitude / (1 - TurnPoint3) / durat * (durat - t)) * ...
            %     sin(fre * t + PHI);
        end
        ind = ind + 1;
    end 
end